iSSF Model Plotting
Here we are plotting the results of integrated step selection function (iSSF) models fitted to water buffalo in Arnhem Land. These models contain the following continuous predictors:
- Monthly Normalised Difference Vegetation Index (NDVI)
- Monthly Normalised Difference Vegetation Index ^ 2 (NDVI^2)
- Monthly Normalised Burn Ratio with monthly lag (either 0, 1, 2 or 3 months of lag) (NBR lag 0,1,2,3)
- Monthly Distance to water (water presence inferred from monthly Normalised Difference Water Index (NDWI) layer)
- Movement covariates to correct for bias in the iSSF model fitting, and can be used to update movement behaviour
- Step length
- Natural logarithm of the step length
- Cosine of the turning angles
Load packages
Import data
Rows: 1062908 Columns: 58
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (4): pheno_start, pheno_end, season, day_period
dbl (51): ...1, id, burst_, x1_, x2_, y1_, y2_, sl_, ta_, dt_, step_id_, sr...
lgl (1): case_
dttm (2): t1_, t2_
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Get scaling factors for the NDVI and NDVI^2 covariates
Code
ssf_data <- ssf_data %>% mutate(srtm_start_scaled = scale(srtm_start),
srtm_end_scaled = scale(srtm_end),
water_dist_start_scaled = scale(water_dist_start),
water_dist_end_scaled = scale(water_dist_end),
ndvi_start_scaled = scale(ndvi_start),
ndvi_end_scaled = scale(ndvi_end),
nbr_start_scaled = scale(nbr_start),
nbr_end_scaled = scale(nbr_end),
sl_scaled = scale(sl_),
log_sl_scaled = scale(log_sl_),
cos_ta_scaled = scale(cos_ta_)
)
# get scaling factors
sl_scale_sd <- attr(ssf_data$sl_scaled, "scaled:scale")
log_sl_scale_sd <- attr(ssf_data$log_sl_scaled, "scaled:scale")
cos_ta_scale_sd <- attr(ssf_data$cos_ta_scaled, "scaled:scale")
ndvi_scale_sd <- attr(ssf_data$ndvi_start_scaled, "scaled:scale")
ndvi_scale_sd <- 0.1262109 # from line of code aboveCheck the NBR layers and NAFI mask
Import layers
Code
nbr <- rast("nbr_raster_amt_cropped_epsg32753.tif")
terra::time(nbr) <- as.POSIXct(lubridate::ymd("2018-01-01") + months(0:23))
# the NAFI layers - monthly fires indicated by the value
NAFI2018_cropped_resampled <- rast("NAFI2018_cropped_resampled_raster_epsg32753.tif")
NAFI2019_cropped_resampled <- rast("NAFI2019_cropped_resampled_raster_epsg32753.tif")
nbr2018_NAFImasked <- rast("nbr2018_NAFImask_epsg32753.tif")
nbr2019_NAFImasked <- rast("nbr2019_NAFImask_epsg32753.tif")
nbr_NAFImask <- c(nbr2018_NAFImasked, nbr2019_NAFImasked)Plot the NAFI layers
Plot the NBR layers and NAFI masked layers
Create a colour scheme
Load fitted model objects
These models were fitted with the glmmTMB package, and have a lag to the NBR covariate of 0, 1, 2 or 3 months.
Code
dry_lag0 <- readRDS("fitted_models/buffalo_ssf10rs_dry_nbrNAFI_lag0_2024-06-13.rds")
dry_lag1 <- readRDS("fitted_models/buffalo_ssf10rs_dry_nbr_NAFImask_lag1_2024-06-13.rds")
dry_lag2 <- readRDS("fitted_models/buffalo_ssf10rs_dry_nbr_NAFImask_lag2_2024-06-13.rds")
dry_lag3 <- readRDS("fitted_models/buffalo_ssf10rs_dry_nbr_NAFImask_lag3_2024-06-13.rds")Code
wet_lag0 <- readRDS("fitted_models/buffalo_ssf10rs_wet_nbrNAFI_lag0_2024-06-13.rds")
wet_lag1 <- readRDS("fitted_models/buffalo_ssf10rs_wet_nbr_NAFImask_lag1_2024-06-13.rds")
wet_lag2 <- readRDS("fitted_models/buffalo_ssf10rs_wet_nbr_NAFImask_lag2_2024-06-13.rds")
wet_lag3 <- readRDS("fitted_models/buffalo_ssf10rs_wet_nbr_NAFImask_lag3_2024-06-13.rds")Check model outputs
Dry season - NBR mask
Family: poisson ( log )
Formula:
case_ ~ -1 + ndvi_end_scaled + I(ndvi_end_scaled^2) + nbr_NAFImask_end_scaled +
water_dist_end_scaled + sl_scaled + log_sl_scaled + cos_ta_scaled +
(1 | step_id) + (0 + ndvi_end_scaled | id) + (0 + I(ndvi_end_scaled^2) |
id) + (0 + nbr_NAFImask_end_scaled | id) + (0 + water_dist_end_scaled |
id)
Data: ssfdat_dry
AIC BIC logLik deviance df.resid
1158463.1 1158587.8 -579220.6 1158441.1 619424
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
step_id (Intercept) 1.000e+06 1000.0000
id ndvi_end_scaled 4.657e-02 0.2158
id.1 I(ndvi_end_scaled^2) 4.288e-02 0.2071
id.2 nbr_NAFImask_end_scaled 1.101e-02 0.1049
id.3 water_dist_end_scaled 1.944e-01 0.4409
Number of obs: 619435, groups: step_id, 56326; id, 15
Conditional model:
Estimate Std. Error z value Pr(>|z|)
ndvi_end_scaled -0.064210 0.056860 -1.129 0.258790
I(ndvi_end_scaled^2) -0.201018 0.054289 -3.703 0.000213 ***
nbr_NAFImask_end_scaled 0.068109 0.031062 2.193 0.028330 *
water_dist_end_scaled -0.534262 0.119417 -4.474 7.68e-06 ***
sl_scaled 0.023810 0.005625 4.233 2.31e-05 ***
log_sl_scaled -0.006310 0.005485 -1.150 0.249948
cos_ta_scaled -0.015570 0.004440 -3.507 0.000454 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
Family: poisson ( log )
Formula:
case_ ~ -1 + ndvi_end_scaled + I(ndvi_end_scaled^2) + nbr_NAFImask_lag_1_end_scaled +
water_dist_end_scaled + sl_scaled + log_sl_scaled + cos_ta_scaled +
(1 | step_id) + (0 + ndvi_end_scaled | id) + (0 + I(ndvi_end_scaled^2) |
id) + (0 + nbr_NAFImask_lag_1_end_scaled | id) + (0 + water_dist_end_scaled |
id)
Data: ssfdat_dry
AIC BIC logLik deviance df.resid
1158413.3 1158538.0 -579195.6 1158391.3 619324
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
step_id (Intercept) 1.000e+06 1.000e+03
id ndvi_end_scaled 4.862e-02 2.205e-01
id.1 I(ndvi_end_scaled^2) 4.140e-02 2.035e-01
id.2 nbr_NAFImask_lag_1_end_scaled 4.793e-03 6.923e-02
id.3 water_dist_end_scaled 1.973e-01 4.442e-01
Number of obs: 619335, groups: step_id, 56326; id, 15
Conditional model:
Estimate Std. Error z value Pr(>|z|)
ndvi_end_scaled -0.059419 0.058059 -1.023 0.306106
I(ndvi_end_scaled^2) -0.196315 0.053372 -3.678 0.000235 ***
nbr_NAFImask_lag_1_end_scaled 0.049287 0.021860 2.255 0.024153 *
water_dist_end_scaled -0.539565 0.120229 -4.488 7.20e-06 ***
sl_scaled 0.024931 0.005623 4.434 9.26e-06 ***
log_sl_scaled -0.007381 0.005480 -1.347 0.178029
cos_ta_scaled -0.015518 0.004440 -3.495 0.000474 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
Family: poisson ( log )
Formula:
case_ ~ -1 + ndvi_end_scaled + I(ndvi_end_scaled^2) + nbr_NAFImask_lag_2_end_scaled +
water_dist_end_scaled + sl_scaled + log_sl_scaled + cos_ta_scaled +
(1 | step_id) + (0 + ndvi_end_scaled | id) + (0 + I(ndvi_end_scaled^2) |
id) + (0 + nbr_NAFImask_lag_2_end_scaled | id) + (0 + water_dist_end_scaled |
id)
Data: ssfdat_dry
AIC BIC logLik deviance df.resid
1158327.5 1158452.2 -579152.8 1158305.5 619211
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
step_id (Intercept) 1.000e+06 1000.0000
id ndvi_end_scaled 5.088e-02 0.2256
id.1 I(ndvi_end_scaled^2) 4.104e-02 0.2026
id.2 nbr_NAFImask_lag_2_end_scaled 1.196e-02 0.1093
id.3 water_dist_end_scaled 2.011e-01 0.4485
Number of obs: 619222, groups: step_id, 56326; id, 15
Conditional model:
Estimate Std. Error z value Pr(>|z|)
ndvi_end_scaled -0.063999 0.059344 -1.078 0.280840
I(ndvi_end_scaled^2) -0.197331 0.053151 -3.713 0.000205 ***
nbr_NAFImask_lag_2_end_scaled -0.009325 0.031048 -0.300 0.763906
water_dist_end_scaled -0.536862 0.121287 -4.426 9.58e-06 ***
sl_scaled 0.024958 0.005622 4.440 9.01e-06 ***
log_sl_scaled -0.006762 0.005483 -1.233 0.217506
cos_ta_scaled -0.015806 0.004440 -3.560 0.000372 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
Family: poisson ( log )
Formula:
case_ ~ -1 + ndvi_end_scaled + I(ndvi_end_scaled^2) + nbr_NAFImask_lag_3_end_scaled +
water_dist_end_scaled + sl_scaled + log_sl_scaled + cos_ta_scaled +
(1 | step_id) + (0 + ndvi_end_scaled | id) + (0 + I(ndvi_end_scaled^2) |
id) + (0 + nbr_NAFImask_lag_3_end_scaled | id) + (0 + water_dist_end_scaled |
id)
Data: ssfdat_dry
AIC BIC logLik deviance df.resid
1154382.3 1154506.9 -577180.1 1154360.3 617135
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
step_id (Intercept) 1.000e+06 1000.0000
id ndvi_end_scaled 4.670e-02 0.2161
id.1 I(ndvi_end_scaled^2) 4.085e-02 0.2021
id.2 nbr_NAFImask_lag_3_end_scaled 2.016e-02 0.1420
id.3 water_dist_end_scaled 1.892e-01 0.4350
Number of obs: 617146, groups: step_id, 56273; id, 15
Conditional model:
Estimate Std. Error z value Pr(>|z|)
ndvi_end_scaled -0.067141 0.056958 -1.179 0.238489
I(ndvi_end_scaled^2) -0.198354 0.053026 -3.741 0.000184 ***
nbr_NAFImask_lag_3_end_scaled -0.025078 0.039013 -0.643 0.520350
water_dist_end_scaled -0.549929 0.117981 -4.661 3.14e-06 ***
sl_scaled 0.025689 0.005626 4.566 4.96e-06 ***
log_sl_scaled -0.005391 0.005496 -0.981 0.326675
cos_ta_scaled -0.015540 0.004449 -3.493 0.000478 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Wet season - NBR mask
Code
Family: poisson ( log )
Formula:
case_ ~ -1 + ndvi_end_scaled + I(ndvi_end_scaled^2) + nbr_NAFImask_end_scaled +
water_dist_end_scaled + sl_scaled + log_sl_scaled + cos_ta_scaled +
(1 | step_id) + (0 + ndvi_end_scaled | id) + (0 + I(ndvi_end_scaled^2) |
id) + (0 + nbr_NAFImask_end_scaled | id) + (0 + water_dist_end_scaled |
id)
Data: ssfdat_wet
AIC BIC logLik deviance df.resid
780451.5 780571.9 -390214.8 780429.5 416401
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
step_id (Intercept) 1.000e+06 1.000e+03
id ndvi_end_scaled 4.079e-02 2.020e-01
id.1 I(ndvi_end_scaled^2) 3.577e-03 5.981e-02
id.2 nbr_NAFImask_end_scaled 6.461e-03 8.038e-02
id.3 water_dist_end_scaled 3.911e-01 6.254e-01
Number of obs: 416412, groups: step_id, 39336; id, 15
Conditional model:
Estimate Std. Error z value Pr(>|z|)
ndvi_end_scaled -0.110719 0.055136 -2.008 0.04463 *
I(ndvi_end_scaled^2) -0.069847 0.016921 -4.128 3.66e-05 ***
nbr_NAFImask_end_scaled 0.032391 0.023856 1.358 0.17455
water_dist_end_scaled -0.326782 0.175689 -1.860 0.06288 .
sl_scaled -0.019121 0.006832 -2.799 0.00513 **
log_sl_scaled 0.089186 0.007292 12.231 < 2e-16 ***
cos_ta_scaled 0.029285 0.005464 5.360 8.32e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
Family: poisson ( log )
Formula:
case_ ~ -1 + ndvi_end_scaled + I(ndvi_end_scaled^2) + nbr_NAFImask_lag_1_end_scaled +
water_dist_end_scaled + sl_scaled + log_sl_scaled + cos_ta_scaled +
(1 | step_id) + (0 + ndvi_end_scaled | id) + (0 + I(ndvi_end_scaled^2) |
id) + (0 + nbr_NAFImask_lag_1_end_scaled | id) + (0 + water_dist_end_scaled |
id)
Data: ssfdat_wet
AIC BIC logLik deviance df.resid
765469.0 765589.1 -382723.5 765447.0 408459
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
step_id (Intercept) 1.000e+06 1.000e+03
id ndvi_end_scaled 3.354e-02 1.831e-01
id.1 I(ndvi_end_scaled^2) 3.271e-03 5.720e-02
id.2 nbr_NAFImask_lag_1_end_scaled 2.537e-03 5.037e-02
id.3 water_dist_end_scaled 3.363e-01 5.799e-01
Number of obs: 408470, groups: step_id, 39172; id, 15
Conditional model:
Estimate Std. Error z value Pr(>|z|)
ndvi_end_scaled -0.105740 0.050571 -2.091 0.03654 *
I(ndvi_end_scaled^2) -0.068177 0.016291 -4.185 2.85e-05 ***
nbr_NAFImask_lag_1_end_scaled 0.083075 0.017153 4.843 1.28e-06 ***
water_dist_end_scaled -0.205804 0.164563 -1.251 0.21108
sl_scaled -0.020515 0.006906 -2.971 0.00297 **
log_sl_scaled 0.090305 0.007383 12.232 < 2e-16 ***
cos_ta_scaled 0.030390 0.005527 5.498 3.84e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
Family: poisson ( log )
Formula:
case_ ~ -1 + ndvi_end_scaled + I(ndvi_end_scaled^2) + nbr_NAFImask_lag_2_end_scaled +
water_dist_end_scaled + sl_scaled + log_sl_scaled + cos_ta_scaled +
(1 | step_id) + (0 + ndvi_end_scaled | id) + (0 + I(ndvi_end_scaled^2) |
id) + (0 + nbr_NAFImask_lag_2_end_scaled | id) + (0 + water_dist_end_scaled |
id)
Data: ssfdat_wet
AIC BIC logLik deviance df.resid
767710.8 767831.0 -383844.4 767688.8 409719
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
step_id (Intercept) 1.000e+06 1.000e+03
id ndvi_end_scaled 3.798e-02 1.949e-01
id.1 I(ndvi_end_scaled^2) 3.488e-03 5.906e-02
id.2 nbr_NAFImask_lag_2_end_scaled 1.282e-02 1.132e-01
id.3 water_dist_end_scaled 2.998e-01 5.475e-01
Number of obs: 409730, groups: step_id, 39270; id, 15
Conditional model:
Estimate Std. Error z value Pr(>|z|)
ndvi_end_scaled -0.103767 0.053428 -1.942 0.052115 .
I(ndvi_end_scaled^2) -0.066154 0.016748 -3.950 7.81e-05 ***
nbr_NAFImask_lag_2_end_scaled 0.121523 0.032968 3.686 0.000228 ***
water_dist_end_scaled -0.177165 0.156458 -1.132 0.257485
sl_scaled -0.017556 0.006867 -2.557 0.010570 *
log_sl_scaled 0.094926 0.007394 12.839 < 2e-16 ***
cos_ta_scaled 0.032153 0.005519 5.825 5.70e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
Family: poisson ( log )
Formula:
case_ ~ -1 + ndvi_end_scaled + I(ndvi_end_scaled^2) + nbr_NAFImask_lag_3_end_scaled +
water_dist_end_scaled + sl_scaled + log_sl_scaled + cos_ta_scaled +
(1 | step_id) + (0 + ndvi_end_scaled | id) + (0 + I(ndvi_end_scaled^2) |
id) + (0 + nbr_NAFImask_lag_3_end_scaled | id) + (0 + water_dist_end_scaled |
id)
Data: ssfdat_wet
AIC BIC logLik deviance df.resid
767616.8 767736.9 -383797.4 767594.8 409701
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
step_id (Intercept) 1.000e+06 1.000e+03
id ndvi_end_scaled 3.876e-02 1.969e-01
id.1 I(ndvi_end_scaled^2) 3.630e-03 6.025e-02
id.2 nbr_NAFImask_lag_3_end_scaled 1.429e-03 3.780e-02
id.3 water_dist_end_scaled 2.790e-01 5.282e-01
Number of obs: 409712, groups: step_id, 39298; id, 15
Conditional model:
Estimate Std. Error z value Pr(>|z|)
ndvi_end_scaled -0.107474 0.053884 -1.995 0.046092 *
I(ndvi_end_scaled^2) -0.065768 0.017046 -3.858 0.000114 ***
nbr_NAFImask_lag_3_end_scaled 0.063826 0.016883 3.780 0.000157 ***
water_dist_end_scaled -0.199873 0.151718 -1.317 0.187705
sl_scaled -0.017101 0.006864 -2.491 0.012723 *
log_sl_scaled 0.092028 0.007384 12.463 < 2e-16 ***
cos_ta_scaled 0.031568 0.005519 5.720 1.07e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Create dataframes to plot the coefficients
Code
coef_df_dry_lag0 <- data.frame("model" = "dry_lag0",
"model_covariate" = names(fixef(dry_lag0)$cond),
"Covariate" = c("NDVI", "NDVI^2", "NBR", "DistanceToWater", "sl", "log(sl)", "cos(ta)"),
"Estimate" = coef(summary(dry_lag0))$cond[, "Estimate"],
"SE" = coef(summary(dry_lag0))$cond[, "Std. Error"]
)
coef_df_dry_lag1 <- data.frame("model" = "dry_lag1",
"model_covariate" = names(fixef(dry_lag1)$cond),
"Covariate" = c("NDVI", "NDVI^2", "NBR", "DistanceToWater", "sl", "log(sl)", "cos(ta)"),
"Estimate" = coef(summary(dry_lag1))$cond[, "Estimate"],
"SE" = coef(summary(dry_lag1))$cond[, "Std. Error"]
)
coef_df_dry_lag2 <- data.frame("model" = "dry_lag2",
"model_covariate" = names(fixef(dry_lag2)$cond),
"Covariate" = c("NDVI", "NDVI^2", "NBR", "DistanceToWater", "sl", "log(sl)", "cos(ta)"),
"Estimate" = coef(summary(dry_lag2))$cond[, "Estimate"],
"SE" = coef(summary(dry_lag2))$cond[, "Std. Error"]
)
coef_df_dry_lag3 <- data.frame("model" = "dry_lag3",
"model_covariate" = names(fixef(dry_lag3)$cond),
"Covariate" = c("NDVI", "NDVI^2", "NBR", "DistanceToWater", "sl", "log(sl)", "cos(ta)"),
"Estimate" = coef(summary(dry_lag3))$cond[, "Estimate"],
"SE" = coef(summary(dry_lag3))$cond[, "Std. Error"]
)
coef_df_dry <- rbind(coef_df_dry_lag0, coef_df_dry_lag1, coef_df_dry_lag2, coef_df_dry_lag3)
coef_df_dry <- coef_df_dry %>%
mutate(
# Calculate confidence intervals
LCI_90 = Estimate - qnorm(1 - (1 - 0.90) / 2) * SE,
UCI_90 = Estimate + qnorm(1 - (1 - 0.90) / 2) * SE,
LCI_95 = Estimate - qnorm(1 - (1 - 0.95) / 2) * SE,
UCI_95 = Estimate + qnorm(1 - (1 - 0.95) / 2) * SE,
LCI_99 = Estimate - qnorm(1 - (1 - 0.99) / 2) * SE,
UCI_99 = Estimate + qnorm(1 - (1 - 0.99) / 2) * SE,
# Compute p-values
pvalue = 2 * pnorm(-abs(Estimate) / SE)
)
# Add stars indicating the significance
coef_df_dry$Significance <- sapply(1:nrow(coef_df_dry), function(x){
if (coef_df_dry$pvalue[x] <= 0.001){
return("***")
} else if (coef_df_dry$pvalue[x] <= 0.01) {
return("**")
} else if (coef_df_dry$pvalue[x] <= 0.05) {
return("*")
}
})
# Remove the intercept term
coef_df_dry <- coef_df_dry %>% filter(Covariate != "Intercept")
# Add a column indicating the preference
coef_df_dry$Preference <- ifelse(coef_df_dry$Estimate > 0, "Preferred", "Avoided")
coef_df_dry$Preference <- factor(coef_df_dry$Preference, levels = c("Avoided", "Preferred"))Prepare data frame for plotting
Specify the order in which the coefficients should be plotted
Prepare dataset for plotting confidence intervals
Code
# Reverse the order of the 'model' factor levels
coef_df_dry$model <- factor(coef_df_dry$model, levels = rev(levels(factor(coef_df_dry$model))))
coef_df_dry2 <- coef_df_dry %>%
dplyr::select(model, Covariate, Estimate, Preference, LCI_90:UCI_99) %>%
gather(key = confidence_level, value = value, LCI_90:UCI_99) %>%
separate(col = confidence_level, into = c("Type", "Level"), sep = "_") %>%
spread(key = Type, value = value) %>%
mutate(Level = paste0(Level, "%"))Prepare plot with covariates on the y-axis and the estimates on the x-axis
Code
dry_nbr_lag_plot <- ggplot(data = coef_df_dry,
aes(y = Covariate, x = Estimate, col = factor(Preference), group = factor(model))) +
geom_vline(xintercept = 0, color = "black", lty = 2, lwd = 0.3) +
annotate(geom = "segment",
x = -0.85, xend = 0.25,
y = 3.5, yend = 3.5,
colour = "gray80", lty = 1, lwd = 0.3) +
geom_point(shape = 3, size = 2.5, position = position_dodge(width = 0.5)) +
geom_errorbarh(data = coef_df_dry2,
aes(xmin = LCI, xmax = UCI, linewidth = factor(Level)),
height = 0, alpha = 0.5, position = position_dodge(width = 0.5)) +
geom_text(aes(label = Significance, hjust = 0.5, vjust = 0),
show.legend = F, position = position_dodge(width = 0.5)) +
scale_y_discrete(limits = rev(order)) +
theme_classic() +
scale_x_continuous(limits = c(-0.9, 0.25)) +
coord_capped_cart(left = "both", bottom = "both") +
labs(x = expression(beta*"-Estimate")) +
scale_color_manual(name = "Preference",
values = c("#FF8123", "#9B4200")) +
scale_linewidth_manual(name = "Confidence Level",
values = c(2, 1, 0.3)) +
theme(legend.position = "bottom",
legend.margin = margin(0, 50, 0, -20),
legend.box.margin = margin(-5, -10, -5, -10),
legend.text = element_text(face = 3),
legend.title = element_text(face = 3)) +
guides(colour = guide_legend(title.position = "top", title.hjust = 0.5),
linewidth = guide_legend(title.position = "top", title.hjust = 0.5,
override.aes = list(colour = "#FF8123")))
dry_nbr_lag_plotWarning: Removed 10 rows containing missing values or values outside the scale range
(`geom_text()`).
Most ‘informative’ model
Code
dry_nbr_lag_plot <- ggplot(data = coef_df_dry %>% filter(model == "dry_lag0"),
aes(y = Covariate, x = Estimate, col = factor(Preference), group = factor(model))) +
geom_vline(xintercept = 0, color = "black", lty = 2, lwd = 0.3) +
annotate(geom = "segment",
x = -0.85, xend = 0.25,
y = 3.5, yend = 3.5,
colour = "gray80", lty = 1, lwd = 0.3) +
geom_point(shape = 3, size = 2.5, position = position_dodge(width = 0.5)) +
geom_errorbarh(data = coef_df_dry2 %>% filter(model == "dry_lag0"),
aes(xmin = LCI, xmax = UCI, linewidth = factor(Level)),
height = 0, alpha = 0.5, position = position_dodge(width = 0.5)) +
geom_text(aes(label = Significance, hjust = 0.5, vjust = 0),
show.legend = F, position = position_dodge(width = 0.5)) +
scale_y_discrete(limits = rev(order)) +
theme_classic() +
scale_x_continuous(limits = c(-0.9, 0.25)) +
coord_capped_cart(left = "both", bottom = "both") +
labs(x = expression(beta*"-Estimate")) +
scale_color_manual(name = "Preference",
values = c("#FF8123", "#9B4200")) +
scale_linewidth_manual(name = "Confidence Level",
values = c(2, 1, 0.3)) +
theme(legend.position = "bottom",
legend.margin = margin(0, 50, 0, -20),
legend.box.margin = margin(-5, -10, -5, -10),
legend.text = element_text(face = 3),
legend.title = element_text(face = 3)) +
guides(colour = guide_legend(title.position = "top", title.hjust = 0.5),
linewidth = guide_legend(title.position = "top", title.hjust = 0.5,
override.aes = list(colour = "#FF8123")))
dry_nbr_lag_plotWarning: Removed 2 rows containing missing values or values outside the scale range
(`geom_text()`).
NDVI quadratic curve
Pull out the coefficients from the fitted model.
Code
best_nbr_model <- coef_df_dry %>% filter(model == "dry_lag0")
ndvi_param <- best_nbr_model %>% filter(Covariate == "NDVI") %>% pull(Estimate) / ndvi_scale_sd
ndvi_param_se <- best_nbr_model %>% filter(Covariate == "NDVI") %>% pull(SE) / ndvi_scale_sd
ndvi2_param <- best_nbr_model %>% filter(Covariate == "NDVI^2") %>% pull(Estimate) / ndvi_scale_sd
ndvi2_param_se <- best_nbr_model %>% filter(Covariate == "NDVI^2") %>% pull(SE) / ndvi_scale_sdCreate a dataframe of x and y values for the quadratic curve.
Plot the quadratic curve
Code
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
Calculate confidence intervals for the curve.
Here we simulate the beta values for NDVI and NDVI^2, and calculate the RSS for each combination of beta values and NDVI values. We then calculate the 95% confidence intervals for the RSS values.
Code
# Number of simulations
n_sim <- 10000
# Simulate beta values
ndvi_sim <- rnorm(n_sim, mean = ndvi_param, sd = ndvi_param_se)
# hist(ndvi_sim)
ndvi2_sim <- rnorm(n_sim, mean = ndvi2_param, sd = ndvi2_param_se)
# hist(ndvi2_sim)
# Define x values
ndvi_values <- seq(-0.2, 0.8, length.out = 100)
# Initialize matrix to store RSS simulations
RSS_sim_matrix <- matrix(NA, nrow = n_sim, ncol = length(ndvi_values))
# Calculate RSS for each simulated beta and x
for (i in 1:length(ndvi_values)) {
RSS_sim_matrix[, i] <- exp(ndvi_sim * ndvi_values[i] + ndvi2_sim * ndvi_values[i]^2)
}
# Calculate confidence intervals for each x
RSS_CI_lower <- apply(RSS_sim_matrix, 2, quantile, probs = 0.025)
RSS_CI_upper <- apply(RSS_sim_matrix, 2, quantile, probs = 0.975)
# Convert the matrix to a data frame
RSS_sim_df <- as.data.frame(RSS_sim_matrix)
# Add column names corresponding to x_values
colnames(RSS_sim_df) <- paste0("x_", ndvi_values)
# Add a simulation ID column
RSS_sim_df$Simulation <- 1:n_sim
# Reshape the data frame from wide to long format
RSS_long_df <- pivot_longer(
RSS_sim_df,
cols = starts_with("x_"),
names_to = "x",
values_to = "RSS_sim"
)
# Remove the "x_" prefix and convert x to numeric
RSS_long_df$x <- as.numeric(sub("x_", "", RSS_long_df$x))
# Calculate mean and confidence intervals for each x
RSS_summary <- RSS_long_df %>%
group_by(x) %>%
summarise(
RSS_mean = mean(RSS_sim),
Lower_CI = quantile(RSS_sim, probs = 0.025),
Lower_CI_50 = quantile(RSS_sim, probs = 0.25),
Upper_CI = quantile(RSS_sim, probs = 0.975),
Upper_CI_50 = quantile(RSS_sim, probs = 0.75)
)
# Plot mean RSS with confidence intervals
ggplot(RSS_summary, aes(x = x, y = RSS_mean)) +
geom_hline(yintercept = 1, color = "black", lty = "dashed", lwd = 0.3) +
geom_ribbon(aes(ymin = Lower_CI, ymax = Upper_CI), fill = "#FF8123", alpha = 0.15) +
geom_ribbon(aes(ymin = Lower_CI_50, ymax = Upper_CI_50), fill = "#FF8123", alpha = 0.15) +
geom_line(color = "#FF8123") +
labs(
# title = "Relative Selection Strength with 50% and 95% Confidence Intervals",
x = "NDVI",
y = "Relative Selection Strength (RSS)"
) +
theme_bw()Code
coef_df_wet_lag0 <- data.frame("model" = "wet_lag0",
"model_covariate" = names(fixef(wet_lag0)$cond),
"Covariate" = c("NDVI", "NDVI^2", "NBR", "DistanceToWater", "sl", "log(sl)", "cos(ta)"),
"Estimate" = coef(summary(wet_lag0))$cond[, "Estimate"],
"SE" = coef(summary(wet_lag0))$cond[, "Std. Error"]
)
coef_df_wet_lag1 <- data.frame("model" = "wet_lag1",
"model_covariate" = names(fixef(wet_lag1)$cond),
"Covariate" = c("NDVI", "NDVI^2", "NBR", "DistanceToWater", "sl", "log(sl)", "cos(ta)"),
"Estimate" = coef(summary(wet_lag1))$cond[, "Estimate"],
"SE" = coef(summary(wet_lag1))$cond[, "Std. Error"]
)
coef_df_wet_lag2 <- data.frame("model" = "wet_lag2",
"model_covariate" = names(fixef(wet_lag2)$cond),
"Covariate" = c("NDVI", "NDVI^2", "NBR", "DistanceToWater", "sl", "log(sl)", "cos(ta)"),
"Estimate" = coef(summary(wet_lag2))$cond[, "Estimate"],
"SE" = coef(summary(wet_lag2))$cond[, "Std. Error"]
)
coef_df_wet_lag3 <- data.frame("model" = "wet_lag3",
"model_covariate" = names(fixef(wet_lag3)$cond),
"Covariate" = c("NDVI", "NDVI^2", "NBR", "DistanceToWater", "sl", "log(sl)", "cos(ta)"),
"Estimate" = coef(summary(wet_lag3))$cond[, "Estimate"],
"SE" = coef(summary(wet_lag3))$cond[, "Std. Error"]
)
coef_df_wet <- rbind(coef_df_wet_lag0, coef_df_wet_lag1, coef_df_wet_lag2, coef_df_wet_lag3)
coef_df_wet <- coef_df_wet %>%
mutate(
# Calculate confidence intervals
LCI_90 = Estimate - qnorm(1 - (1 - 0.90) / 2) * SE,
UCI_90 = Estimate + qnorm(1 - (1 - 0.90) / 2) * SE,
LCI_95 = Estimate - qnorm(1 - (1 - 0.95) / 2) * SE,
UCI_95 = Estimate + qnorm(1 - (1 - 0.95) / 2) * SE,
LCI_99 = Estimate - qnorm(1 - (1 - 0.99) / 2) * SE,
UCI_99 = Estimate + qnorm(1 - (1 - 0.99) / 2) * SE,
# Compute p-values
pvalue = 2 * pnorm(-abs(Estimate) / SE)
)
# Add stars indicating the significance
coef_df_wet$Significance <- sapply(1:nrow(coef_df_wet), function(x){
if (coef_df_wet$pvalue[x] <= 0.001){
return("***")
} else if (coef_df_wet$pvalue[x] <= 0.01) {
return("**")
} else if (coef_df_wet$pvalue[x] <= 0.05) {
return("*")
}
})
# Remove the intercept term
coef_df_wet <- coef_df_wet %>% filter(Covariate != "Intercept")
# Add a column indicating the preference
coef_df_wet$Preference <- ifelse(coef_df_wet$Estimate > 0, "Preferred", "Avoided")
coef_df_wet$Preference <- factor(coef_df_wet$Preference, levels = c("Avoided", "Preferred"))Prepare confidence interval dataframe
Code
# Reverse the order of the 'model' factor levels
coef_df_wet$model <- factor(coef_df_wet$model, levels = rev(levels(factor(coef_df_wet$model))))
coef_df_wet2 <- coef_df_wet %>%
dplyr::select(model, Covariate, Estimate, Preference, LCI_90:UCI_99) %>%
gather(key = confidence_level, value = value, LCI_90:UCI_99) %>%
separate(col = confidence_level, into = c("Type", "Level"), sep = "_") %>%
spread(key = Type, value = value) %>%
mutate(Level = paste0(Level, "%"))Prepare plot with covariates on the y-axis and the estimates on the x-axis
Code
wet_nbr_lag_plot <- ggplot(data = coef_df_wet,
aes(y = Covariate, x = Estimate, col = factor(Preference), group = factor(model))) +
geom_vline(xintercept = 0, color = "black", lty = 2, lwd = 0.3) +
annotate(geom = "segment",
x = -0.9, xend = 0.25,
y = 3.5, yend = 3.5,
colour = "gray80", lty = 1, lwd = 0.3) +
geom_point(shape = 3, size = 2.5, position = position_dodge(width = 0.5)) +
geom_errorbarh(data = coef_df_wet2,
aes(xmin = LCI, xmax = UCI, linewidth = factor(Level)),
height = 0, alpha = 0.5, position = position_dodge(width = 0.5)) +
geom_text(aes(label = Significance, hjust = 0.5, vjust = 0),
show.legend = F, position = position_dodge(width = 0.5)) +
scale_y_discrete(limits = rev(order)) +
theme_classic() +
scale_x_continuous(limits = c(-0.9, 0.25)) +
coord_capped_cart(left = "both", bottom = "both") +
labs(x = expression(beta*"-Estimate")) +
scale_color_manual(name = "Preference",
values = c("#238DFF", "#004691")) +
scale_linewidth_manual(name = "Confidence Level",
values = c(2, 1, 0.3)) +
theme(legend.position = "bottom",
legend.margin = margin(0, 50, 0, -20),
legend.box.margin = margin(-5, -10, -5, -10),
legend.text = element_text(face = 3),
legend.title = element_text(face = 3)) +
guides(colour = guide_legend(title.position = "top", title.hjust = 0.5),
linewidth = guide_legend(title.position = "top", title.hjust = 0.5,
override.aes = list(colour = "#1377E2")))
wet_nbr_lag_plotWarning: Removed 6 rows containing missing values or values outside the scale range
(`geom_text()`).
Most ‘informative’ model
Code
wet_nbr_lag_plot <- ggplot(data = coef_df_wet %>% filter(model == "wet_lag2"),
aes(y = Covariate, x = Estimate, col = factor(Preference), group = factor(model))) +
geom_vline(xintercept = 0, color = "black", lty = 2, lwd = 0.3) +
annotate(geom = "segment",
x = -0.9, xend = 0.25,
y = 3.5, yend = 3.5,
colour = "gray80", lty = 1, lwd = 0.3) +
geom_point(shape = 3, size = 2.5, position = position_dodge(width = 0.5)) +
geom_errorbarh(data = coef_df_wet2 %>% filter(model == "wet_lag2"),
aes(xmin = LCI, xmax = UCI, linewidth = factor(Level)),
height = 0, alpha = 0.5, position = position_dodge(width = 0.5)) +
geom_text(aes(label = Significance, hjust = 0.5, vjust = 0),
show.legend = F, position = position_dodge(width = 0.5)) +
scale_y_discrete(limits = rev(order)) +
theme_classic() +
scale_x_continuous(limits = c(-0.9, 0.25)) +
coord_capped_cart(left = "both", bottom = "both") +
labs(x = expression(beta*"-Estimate")) +
scale_color_manual(name = "Preference",
values = c("#238DFF", "#004691")) +
scale_linewidth_manual(name = "Confidence Level",
values = c(2, 1, 0.3)) +
theme(legend.position = "bottom",
legend.margin = margin(0, 50, 0, -20),
legend.box.margin = margin(-5, -10, -5, -10),
legend.text = element_text(face = 3),
legend.title = element_text(face = 3)) +
guides(colour = guide_legend(title.position = "top", title.hjust = 0.5),
linewidth = guide_legend(title.position = "top", title.hjust = 0.5,
override.aes = list(colour = "#1377E2")))
wet_nbr_lag_plotWarning: Removed 2 rows containing missing values or values outside the scale range
(`geom_text()`).
NDVI quadratic curve
Pull out the coefficients from the fitted model.
Code
best_nbr_model <- coef_df_wet %>% filter(model == "wet_lag0")
ndvi_param <- best_nbr_model %>% filter(Covariate == "NDVI") %>% pull(Estimate) / ndvi_scale_sd
ndvi_param_se <- best_nbr_model %>% filter(Covariate == "NDVI") %>% pull(SE) / ndvi_scale_sd
ndvi2_param <- best_nbr_model %>% filter(Covariate == "NDVI^2") %>% pull(Estimate) / ndvi_scale_sd
ndvi2_param_se <- best_nbr_model %>% filter(Covariate == "NDVI^2") %>% pull(SE) / ndvi_scale_sdCreate a dataframe of x and y values for the quadratic curve.
Plot the quadratic curve.
Code
Calculate confidence intervals for the curve.
Here we simulate the beta values for NDVI and NDVI^2, and calculate the RSS for each combination of beta values and NDVI values. We then calculate the 95% confidence intervals for the RSS values.
Code
# Number of simulations
n_sim <- 10000
# Simulate beta values
ndvi_sim <- rnorm(n_sim, mean = ndvi_param, sd = ndvi_param_se)
# hist(ndvi_sim)
ndvi2_sim <- rnorm(n_sim, mean = ndvi2_param, sd = ndvi2_param_se)
# hist(ndvi2_sim)
# Define x values
ndvi_values <- seq(-0.2, 0.8, length.out = 100)
# Initialize matrix to store RSS simulations
RSS_sim_matrix <- matrix(NA, nrow = n_sim, ncol = length(ndvi_values))
# Calculate RSS for each simulated beta and x
for (i in 1:length(ndvi_values)) {
RSS_sim_matrix[, i] <- exp(ndvi_sim * ndvi_values[i] + ndvi2_sim * ndvi_values[i]^2)
}
# Calculate confidence intervals for each x
RSS_CI_lower <- apply(RSS_sim_matrix, 2, quantile, probs = 0.025)
RSS_CI_upper <- apply(RSS_sim_matrix, 2, quantile, probs = 0.975)
# Convert the matrix to a data frame
RSS_sim_df <- as.data.frame(RSS_sim_matrix)
# Add column names corresponding to x_values
colnames(RSS_sim_df) <- paste0("x_", ndvi_values)
# Add a simulation ID column
RSS_sim_df$Simulation <- 1:n_sim
# Reshape the data frame from wide to long format
RSS_long_df <- pivot_longer(
RSS_sim_df,
cols = starts_with("x_"),
names_to = "x",
values_to = "RSS_sim"
)
# Remove the "x_" prefix and convert x to numeric
RSS_long_df$x <- as.numeric(sub("x_", "", RSS_long_df$x))
# Calculate mean and confidence intervals for each x
RSS_summary <- RSS_long_df %>%
group_by(x) %>%
summarise(
RSS_mean = mean(RSS_sim),
Lower_CI = quantile(RSS_sim, probs = 0.025),
Lower_CI_50 = quantile(RSS_sim, probs = 0.25),
Upper_CI = quantile(RSS_sim, probs = 0.975),
Upper_CI_50 = quantile(RSS_sim, probs = 0.75)
)
# Plot mean RSS with confidence intervals
ggplot(RSS_summary, aes(x = x, y = RSS_mean)) +
geom_hline(yintercept = 1, color = "black", lty = "dashed", lwd = 0.3) +
geom_ribbon(aes(ymin = Lower_CI, ymax = Upper_CI), fill = "#1377E2", alpha = 0.15) +
geom_ribbon(aes(ymin = Lower_CI_50, ymax = Upper_CI_50), fill = "#1377E2", alpha = 0.15) +
geom_line(color = "#1377E2") +
labs(
# title = "Relative Selection Strength with 50% and 95% Confidence Intervals",
x = "NDVI",
y = "Relative Selection Strength (RSS)"
) +
theme_bw()Checking used vs avilable distributions
NDVI
Code
# bin the ssf_data$nbr_end values
# ssf_data$ndvi_end_bin <- cut(ssf_data$ndvi_end, breaks = c(-1, -0.5, -0.25, -0.1, 0.1, 0.25, 0.5, 1))
ssf_data$ndvi_end_bin <- cut(ssf_data$ndvi_end, breaks = seq(-1,1,0.1))
ssf_data_plot <- ssf_data %>% filter(ndvi_end > -0.2 & ndvi_end < 0.8)
# inspect the data for dry season
ssf_data_plot %>%
filter(season == "dry") %>%
drop_na() %>%
group_by(case_, ndvi_end_bin) %>%
summarize(n = n()) %>%
mutate(prop = n / sum(n),
label = paste0(round(prop * 100, 1), "%")) %>%
ggplot(aes(ndvi_end_bin, prop, fill = case_, group=case_,label = label)) +
geom_col(position = position_dodge2()) +
geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
labs(x = "NDVI", y = "Proportion", fill = "Case")+
ggtitle("Dry season")+
scale_fill_brewer(palette = "Paired", name="case_",
breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
theme_classic() +
theme(axis.text.x = element_text(angle = 20, hjust = 1))`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
Code
# inspect the data for wet season
ssf_data_plot %>%
filter(season == "wet") %>%
drop_na() %>%
group_by(case_, ndvi_end_bin) %>%
summarize(n = n()) %>%
mutate(prop = n / sum(n),
label = paste0(round(prop * 100, 1), "%")) %>%
ggplot(aes(ndvi_end_bin, prop, fill = case_, group=case_,label = label)) +
geom_col(position = position_dodge2()) +
geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
labs(x = "NDVI", y = "Proportion", fill = "Case")+
ggtitle("Wet season")+
scale_fill_brewer(palette = "Paired", name="case_",
breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
theme_classic() +
theme(axis.text.x = element_text(angle = 20, hjust = 1))`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
NBR
Code
# bin the ssf_data$nbr_end values
ssf_data$nbr_lag_0_end_bin <- cut(ssf_data$nbr_lag_0_end, breaks = c(-1, -0.5, -0.25, -0.1, 0, 0.1, 0.25, 0.5, 1))
# inspect the data for dry season
ssf_data %>%
filter(season == "dry") %>%
drop_na() %>%
group_by(case_, nbr_lag_0_end_bin) %>%
summarize(n = n()) %>%
mutate(prop = n / sum(n),
label = paste0(round(prop * 100, 1), "%")) %>%
ggplot(aes(nbr_lag_0_end_bin, prop, fill = case_, group=case_,label = label)) +
geom_col(position = position_dodge2()) +
geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
labs(x = "Land use class", y = "Proportion", fill = "Case")+
ggtitle("Dry season - no NBR lag")+
scale_fill_brewer(palette = "Paired", name="case_",
breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
theme_classic() +
theme(axis.text.x = element_text(angle = 20, hjust = 1))`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
Code
# inspect the data for wet season
ssf_data %>%
filter(season == "wet") %>%
drop_na() %>%
group_by(case_, nbr_lag_0_end_bin) %>%
summarize(n = n()) %>%
mutate(prop = n / sum(n),
label = paste0(round(prop * 100, 1), "%")) %>%
ggplot(aes(nbr_lag_0_end_bin, prop, fill = case_, group=case_,label = label)) +
geom_col(position = position_dodge2()) +
geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
labs(x = "Land use class", y = "Proportion", fill = "Case")+
ggtitle("Wet season - no NBR lag")+
scale_fill_brewer(palette = "Paired", name="case_",
breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
theme_classic() +
theme(axis.text.x = element_text(angle = 20, hjust = 1))`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
Code
# bin the ssf_data$nbr_end values
ssf_data$nbr_lag_1_end_bin <- cut(ssf_data$nbr_lag_1_end, breaks = c(-1, -0.5, -0.25, -0.1, 0, 0.1, 0.25, 0.5, 1))
# inspect the data for dry season
ssf_data %>%
filter(season == "dry") %>%
drop_na() %>%
group_by(case_, nbr_lag_1_end_bin) %>%
summarize(n = n()) %>%
mutate(prop = n / sum(n),
label = paste0(round(prop * 100, 1), "%")) %>%
ggplot(aes(nbr_lag_1_end_bin, prop, fill = case_, group=case_,label = label)) +
geom_col(position = position_dodge2()) +
geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
labs(x = "Land use class", y = "Proportion", fill = "Case")+
ggtitle("Dry season - 1 month NBR lag")+
scale_fill_brewer(palette = "Paired", name="case_",
breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
theme_classic() +
theme(axis.text.x = element_text(angle = 20, hjust = 1))`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
Code
# inspect the data for wet season
ssf_data %>%
filter(season == "wet") %>%
drop_na() %>%
group_by(case_, nbr_lag_1_end_bin) %>%
summarize(n = n()) %>%
mutate(prop = n / sum(n),
label = paste0(round(prop * 100, 1), "%")) %>%
ggplot(aes(nbr_lag_1_end_bin, prop, fill = case_, group=case_,label = label)) +
geom_col(position = position_dodge2()) +
geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
labs(x = "Land use class", y = "Proportion", fill = "Case")+
ggtitle("Wet season - 1 month NBR lag")+
scale_fill_brewer(palette = "Paired", name="case_",
breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
theme_classic() +
theme(axis.text.x = element_text(angle = 20, hjust = 1))`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
Code
# bin the ssf_data$nbr_end values
ssf_data$nbr_lag_2_end_bin <- cut(ssf_data$nbr_lag_2_end, breaks = c(-1, -0.5, -0.25, -0.1, 0, 0.1, 0.25, 0.5, 1))
# inspect the data for dry season
ssf_data %>%
filter(season == "dry") %>%
drop_na() %>%
group_by(case_, nbr_lag_2_end_bin) %>%
summarize(n = n()) %>%
mutate(prop = n / sum(n),
label = paste0(round(prop * 100, 1), "%")) %>%
ggplot(aes(nbr_lag_2_end_bin, prop, fill = case_, group=case_,label = label)) +
geom_col(position = position_dodge2()) +
geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
labs(x = "Land use class", y = "Proportion", fill = "Case")+
ggtitle("Dry season - 2 month NBR lag")+
scale_fill_brewer(palette = "Paired", name="case_",
breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
theme_classic() +
theme(axis.text.x = element_text(angle = 20, hjust = 1))`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
Code
# inspect the data for wet season
ssf_data %>%
filter(season == "wet") %>%
drop_na() %>%
group_by(case_, nbr_lag_2_end_bin) %>%
summarize(n = n()) %>%
mutate(prop = n / sum(n),
label = paste0(round(prop * 100, 1), "%")) %>%
ggplot(aes(nbr_lag_2_end_bin, prop, fill = case_, group=case_,label = label)) +
geom_col(position = position_dodge2()) +
geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
labs(x = "Land use class", y = "Proportion", fill = "Case")+
ggtitle("Wet season - 2 month NBR lag")+
scale_fill_brewer(palette = "Paired", name="case_",
breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
theme_classic() +
theme(axis.text.x = element_text(angle = 20, hjust = 1))`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
Code
# bin the ssf_data$nbr_end values
ssf_data$nbr_lag_3_end_bin <- cut(ssf_data$nbr_lag_3_end, breaks = c(-1, -0.5, -0.25, -0.1, 0, 0.1, 0.25, 0.5, 1))
# inspect the data for dry season
ssf_data %>%
filter(season == "dry") %>%
drop_na() %>%
group_by(case_, nbr_lag_3_end_bin) %>%
summarize(n = n()) %>%
mutate(prop = n / sum(n),
label = paste0(round(prop * 100, 1), "%")) %>%
ggplot(aes(nbr_lag_3_end_bin, prop, fill = case_, group=case_,label = label)) +
geom_col(position = position_dodge2()) +
geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
labs(x = "Land use class", y = "Proportion", fill = "Case")+
ggtitle("Dry season - 3 month NBR lag")+
scale_fill_brewer(palette = "Paired", name="case_",
breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
theme_classic() +
theme(axis.text.x = element_text(angle = 20, hjust = 1))`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
Code
# inspect the data for wet season
ssf_data %>%
filter(season == "wet") %>%
drop_na() %>%
group_by(case_, nbr_lag_3_end_bin) %>%
summarize(n = n()) %>%
mutate(prop = n / sum(n),
label = paste0(round(prop * 100, 1), "%")) %>%
ggplot(aes(nbr_lag_3_end_bin, prop, fill = case_, group=case_,label = label)) +
geom_col(position = position_dodge2()) +
geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
labs(x = "Land use class", y = "Proportion", fill = "Case")+
ggtitle("Wet season - 3 month NBR lag")+
scale_fill_brewer(palette = "Paired", name="case_",
breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
theme_classic() +
theme(axis.text.x = element_text(angle = 20, hjust = 1))`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
R version 4.4.1 (2024-06-14 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 10 x64 (build 19045)
Matrix products: default
locale:
[1] LC_COLLATE=English_Australia.utf8 LC_CTYPE=English_Australia.utf8
[3] LC_MONETARY=English_Australia.utf8 LC_NUMERIC=C
[5] LC_TIME=English_Australia.utf8
time zone: Australia/Brisbane
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] lemon_0.4.9 RColorBrewer_1.1-3 terra_1.7-78 glmmTMB_1.1.10
[5] lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4
[9] purrr_1.0.2 readr_2.1.5 tidyr_1.3.1 tibble_3.2.1
[13] ggplot2_3.5.1 tidyverse_2.0.0
loaded via a namespace (and not attached):
[1] gtable_0.3.5 TMB_1.9.15 xfun_0.47
[4] htmlwidgets_1.6.4 lattice_0.22-6 tzdb_0.4.0
[7] numDeriv_2016.8-1.1 vctrs_0.6.5 tools_4.4.1
[10] Rdpack_2.6.1 generics_0.1.3 parallel_4.4.1
[13] fansi_1.0.6 pkgconfig_2.0.3 Matrix_1.7-0
[16] lifecycle_1.0.4 farver_2.1.2 compiler_4.4.1
[19] textshaping_0.4.0 munsell_0.5.1 codetools_0.2-20
[22] htmltools_0.5.8.1 yaml_2.3.10 crayon_1.5.3
[25] pillar_1.9.0 nloptr_2.1.1 MASS_7.3-64
[28] reformulas_0.3.0 boot_1.3-30 nlme_3.1-164
[31] tidyselect_1.2.1 digest_0.6.37 mvtnorm_1.3-3
[34] stringi_1.8.4 labeling_0.4.3 splines_4.4.1
[37] fastmap_1.2.0 grid_4.4.1 colorspace_2.1-1
[40] cli_3.6.3 magrittr_2.0.3 utf8_1.2.4
[43] withr_3.0.1 scales_1.3.0 bit64_4.0.5
[46] timechange_0.3.0 estimability_1.5.1 rmarkdown_2.28
[49] emmeans_1.10.7 bit_4.0.5 lme4_1.1-35.5
[52] gridExtra_2.3 ragg_1.3.3 hms_1.1.3
[55] coda_0.19-4.1 evaluate_1.0.0 knitr_1.48
[58] rbibutils_2.2.16 mgcv_1.9-1 rlang_1.1.4
[61] Rcpp_1.0.13 xtable_1.8-4 glue_1.7.0
[64] vroom_1.6.5 rstudioapi_0.16.0 minqa_1.2.8
[67] jsonlite_1.8.8 plyr_1.8.9 R6_2.5.1
[70] systemfonts_1.1.0